AD-A071 092 
UNCLASSIFIED 


NAVAL SURFACE WEAPONS CENTER DAHLOREN LAB VA f/% a /3 

FIVE M-*yiTHMS FOR THE EARTH TIDE TERNS IN THE FORCE MOOCL OF —ETCfU) 

FEB 79 M OROEBER 

NSWC/DL-TR-391* 


/or | I 

*#71095 

■ 



B 

i 

& 



miXm 

is 



r 






i - v 

s=i 

W- 

i : 







-r- 





'^i_i 










.... 

Ti 














f 






m 



_ 










_ 

EH 8 -79 































© 







UNCLASSIFIED 


V CLASSIFICATION OF THIS P 


pTaEpoRTTr 
NSWC/DL-TR-3912 / 


REPORT DOCUMENTATION PAGE 


GOVT ACCESSION NO. 


RECIPIENT'S CAT Al 


JIVE ALGORITHMS FOR THE EARTH TIDE 

Jerms in the force model of nswc 
COMPUTER PROGRAMS FOR SATELLITE 
&EODESY. ^ _ 




TYPE OF REPORT A PERIOD COVERED 


Final r^X 


SOW. REPORT NUMBER 


S. CONTRACT OR GRANT NUMBERfc; 


Walter/ Grc 


9. PERFORMING ORGANIZATION NAME AND ADORESS 

Naval Surface Weapons Center (K11) 
Dahlgren Laboratory 
Dahlgren, VA 22448 


/ 


PROGRAM ELEMENT. PROJECT. T 
AREA A WORK UNIT NUMBERS 


TZ). 


. CONTROLLING OFFICE NAME AND ADDRESS 


Headquarters 
Defense Mapping Agency 
Washington, DC 20360 


Febnaiy W79 (' 


. NUMBER OF PAGES 


79 


l • r 


15. SECURITY CLASS, (of 

UNCLASSIFIED 

ISn. DECLASSI FI CATION/DOWN GRADING 
SCHEDULE 


IS. DISTRIBUTION STATEMENT (ol t 


Approved for public release; distribution unlimited 


. DISTRIBUTION STATEMENT (at 


. SUPPLEMENTARY NOTES 


19. KEY WORDS rCon 


Satellite geodesy 

Force equation for satellite motion 
Lunar air tide 
Solar air tide 


id Identity by block number) 

Ocean tide 

Tides of the solid earth 
Tidal lag 


*0. ABSTRACT CConl/nu. on rovoroo .M. II nocoxary and Ido nitty by Slock number; 

"^The most recent versions of five computer algorithms which jointly specify the gravitational 
action by which the tidal dislocation of the earth’s masses perturbs satellite motion are presented. 

The first 6f the five>computer routine* formulates the perturbing acceleration which the 
semidiurnal lunar component of the atmospheric tide bulge exerts on the orbit of an artificial earth 
satellite.-^Hie tide bulge is assumed to result from the fact that the earth rotates within the field of 
lunar ma^ attraction, the latter being inhomogeneous across the terrestrial globe. (see back) 




DO , 


1473 7 COITION Of 1 

f S/N 0102-LF.014-6601 


UNCLASSIFIED 


y 9/ 


SECURITY CLASSIFICATION OF THIS PAGE (Who* Dole Knlered) 










UNCLASSIFIED 


SECURITY CLASSIFICATION OF THIS PAGE fWh«n Data Bnlar.dj 


( 20 ) 

7>The second routine augments the first by spelling out the combined effect of solar heating and 
solar mass attractionj The semidiurnal term as well as the diurnal term are considered. 

^The third algorithm accounts for the ocean tide. Only the semidiurnal lunar tide is presently 
considered, but the algorithm sets the pattern for introducing various other spectral components of 
the fluid tide, if desired. The algorithm starts from a table for the global tide amplitudes and phase 
angles. It subsequently approximates the gravitational potential exterior to the earth caused by the 
tidal bulge by that of a system of point masses located on the ocean surface. 

The fourth tide routine is optional; Jti may be used in place of the third algorithm if one wishes 
to base the tide potential on a functionalization of the tidal sea surface. 

Finally, there is a computer routine for the Newtonian attraction caused by the tide bulge of the 
solid earth. Both the lunar and solar tide components are included. The tidal properties are assumed 
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FOREWORD 


Recent progress in instrumentation for data acquisition and computing methods suggests 
that satellite geodesy must concern itself with even such minute effects as the gravitational 
action that the tidal dislocations of the earth’s masses exert upon an artificial satellite. A 
number of computer algorithms suited to this purpose were developed at the Naval Surface 
Weapons Center (NSWC) during the past few years. The relevant physics background as well 
as sketches of the computer routines were published in a sequence of four technical reports 
in 1976 and 1977. Since then, the algorithms were worked out in detail, coded for the 
computer, and subjected to the customary numerical checkout. Also, an alternate procedure 
was developed for the ocean tides, which increased the number of tide, routines from the 
original four to five. The necessary derivations for the additional routine were never 
documented. They are included as Appendix A. The present report is intended to serve as 
a manual of the five algorithms, each being presented in a form ready for immediate 
computer implementation, complete with checkout data. 

The work documented here was done in the Astronautics and Geodesy Division and 
was funded as part of the development of the TERRA computer program which is used to 
determine the earth’s gravity field from terrestrial and satellite observations. Henry E. 
Castro of the Computer Program Division coded those algorithms that involve the 
atmospheric tides and the two versions of the ocean tide. R. Gordon Barker, also of the 
Computer Program Division, coded the solid earth tide. Additionally, these two gentlemen 
invested a great effort while making the required checkout runs on the computer. Raymond 
B. Manrique of the Astronautics and Geodesy Division and the author jointly outlined the 
computer algorithms, with the author bearing the overall responsibility for the finished 
formulations and the sole responsibility for the numerical checkout calculations. 

This report was reviewed by Richard J. Anderle, Head, Astronautics and Geodesy 
Division. 


Released by: 

R. A. NIEMANN, Head 
Strategic Systems Department 
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INTRODUCTION 


Several years ago, while compiling a force model on which to base the trajectory 
integration in the TERRA 1 and CELEST 2 computer programs for satellite geodesy and 
satellite ephemeris computation, the author faced the nearly complete lack of published 
algorithms that would in a suitable form express the gravitational action by which the tidal 
dislocation of the earth’s masses perturbs satellite motion. At the same time, it was realized 
that the earth tides would have to be taken into account if TERRA and CELEST were to 
substantially exceed the accuracy requirements of their forerunners, GEO 3 and ASTRO. 4 
An effort was consequently made to derive the necessary disturbing potentials and 
perturbing accelerations. The result consisted of a number of mathematical models, each 
designed to be realistic and describing one type of perturbation corresponding to one of 
the four earth tides. The details of these models were documented in four technical 
papers. 5,6,7 - 8 Simultaneously, the necessary computer routines were developed and 

submitted for coding. 

The first tide routine formulated the perturbing acceleration which the semidiurnal 

lunar component of the atmospheric tide bulge exerts on the orbit of an artifical satellite. 
The tide bulge was assumed to result from the fact that the earth rotates within the field 

of lunar mass attraction which is inhomogeneous across the terrestrial globe. The routine 

was restricted to the main term of the semidiurnal tide. 

The second tide routine was designed to augment the first by spelling out the 
combined effect on the atmosphere of solar heating and solar mass attraction, plus the 
resulting modification of the earth’s gravitation. The semidiurnal term as well as the diurnal 
term were included. 

The third algorithm was to account for the ocean tide. It started from an existing 

table for the global tide amplitudes and phase angles, and postulated that each tidal height 

value resembles a gravitating point mass that by its presence perturbs the satellite motion. 
Only the semidiurnal lunar tide was considered; however, this algorithm may be regarded as 
setting the pattern for introducing various other spectral components of the fluid tide, if 
desired. 

Further, there was a computer routine for the Newtonian attraction caused by the tide 

bulge of the solid earth. Both the lunar and the solar tide components were included. The 

tidal properties were assumed to be functions of latitude. Accordingly, Love coefficients 
were featured that are zonal harmonic expansions of latitude. An allowance was made for 
tidal lag manifesting itself as a time delay in the response of the tidal mass redistribution 
to the tidal stress field. 
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Several months after completing the coding requests for these four algorithms, an 
opportunity presented itself to devise an interesting alternative to the just mentioned ocean 
tide routine. Namely, a functionalization became available for the tidal ocean surface. This 
consisted of a representation of the tidal ocean surface by an expansion, in surface 
harmonics, that had been least squares fit adjusted to the table of discrete tide amplitudes 
and phases that had been the basis of the ocean tide algorithm. Upon receiving a card 
deck containing the coefficients for this expansion [made available to us through the 
kindness of Dr. C. C. Goad of the National Oceanic and Atmospheric Administration 
(NOAA)l, an optional routine was written, which interprets the functionalized ocean surface 
as a continuous, gravitating, tidal mass layer. From it, the perturbing acceleration acting on 
the satellite, due to the presence of the oceanic tide bulge, is then derived via the Poisson 
integral and subsequent gradient formation. While the physics background and the 
mathematical derivations for the original four tide routines are available in the shape of 
formal reports, 5,6,7,8 the background for this latter optional version of the ocean tide 
algorithm had never before been documented. Thus, to preserve it, it is included in this 
report as Appendix A. 


1 - COMPUTER ROUTINE FOR THE PERTURBING ACCELERATION 
CAUSED BY THE LUNAR AIR TIDE 


INPUT DATA 

R = any reasonable value for the mean radius of the earth, in meters 
G = gravitational constant, in meters, kilograms, and seconds 
suggested value: G = 6.6732E - 11 m 3 kg 1 sec' 2 
A 2 = a constant associated with the physics of the lunar air tide 
suggested value: A 2 = 0.564 kg m' 2 

The parameters J, n, and t* specify the time instant for which the perturbing acceleration 
due to the tide is to be found. Normally, this is the time instant associated with the 
current time line of the orbit integration. 

J = number of the calendar year 

n = number of the day within the year 

t* = mean solar time at Greenwich (GMT, UTC), in seconds 

X]) inertial, Cartesian components of satellite position vector, in meters, 
x 2 |= associated with the just specified time instant (for a precise definition of 
x 3 ) the reference frame, see Appendices B and C of Reference 2) 
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ALGORITHM 


6(ij)=( if (1-1) 

lo i *j 

N = n 6(J, 1975) + (365 + n) 5(J, 1976) + (731 + n) 5(J, 1977) + (1096 + n) 5(J, 1978) 

+ (1461+n)6(J, 1979) +(1826+ n)5(J, 1980) (1-2) 

AT = 5.28E - 04 + (3.56E - 08)N (1-3) 

d = 27392.5 + N + AT + t* (M) 

T = d/36525 (1-5) 

s = 270.434358 + 481267.883141T- 0.001133T 2 + 0.000002T 3- (1-6) 

h = 279.69668 + 36000.768930T + 0.000303T 2 (1-7) 


s and h will result in decimal fractions of degrees. 

Calculate now the earth-fixed, Cartesian satellite coordinates (y 2 , y 2 , and y 3 ) from the 
corresponding inertial coordinates (x p x 2 , and x ? ). 


<ABCD),„,.^x^ (,-8) 


where (ABCD) J n t , is the transformation matrix that manages the transition from the 
inertial frame of the satellite equations of motion (Basic Inertial System, which is 
associated with either 1950.0 or with the beginning of the day UT of trajectory epoch) to 
the earth-fixed reference frame corresponding to the time instant, J, n, t*. This 
transformation is a computer routine that is already part of TERRA/CELEST. It is external 
to the present algorithm. For details, see Appendices B and C of Reference 2. 
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dU _ 
3r 

3r3 ^r5 

- a — 4 - Pj(sin 0) cos 2a + b P 2 (sin 0) cos 2a 

(1-19) 

3U 

ax 

2R 3 2R 5 

- a-jy P 2 (sin 0) sin 2a + b -jy P 2 (sin 0) sin 2a 

( 1 - 20 ) 

au 

30 

3r3 15R 5 

- a -jy- sin 20 cos 2a + b —jy- (7 sin 2 0-4) sin 20 cos 2a 

( 1 - 21 ) 

T = 

+yjy 2 l + y 2 +y 2 

( 1 - 22 ) 

3 = 

a 2 gr 5 -^ 

2 64 

(1-23) 

b = 

Sir 2 a 

A^GR- — — 

2 3072 48 

(1-24) 

P^(x) = 

3(1 - x 2 ) 

(1-25) 

P?(x) = 

15 , , 

— (1 - x 2 )(7x 2 - 1) 

(1-26) 

X = 

sin 0 = 21 
r 

(1-27) 

sin 20 = 

2y 3 Vvi + y 2 

r 2 

(1-28) 

cos 2a = 

cos 2 a - sin 2 a 

(1-29) 

sin 2a = 

2 sin a cos a 

(1-30) 

cos a = 

cos a* cos X - sina* sin X 

(1-31) 

sin a = 

sin a* cos X + cos a* sin X 

(1-32) 

a* = 

t**-p-7?5 

(1-33) 

t** = 

360 ^ 

86400 * 

(1-34) 

v = s 

1 - h 

(1-35) 

cosX = 

y 2 

Vyf +y \ 

(1-36) 
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sin X = 


(1-37) 


V2 

Vyf+yf 


Finally, transform the perturbing acceleration back into the inertial frame 




where (ABCD) T is the transpose of (ABCD). 


(1-38) 


COMPUTER PROGRAM OUTPUT 

The inertial components (T Xj , T X2> and T Xj ) of the perturbing acceleration due to 
the presence of the lunar air tide are output in terms of m sec -2 . If required, convert to 
km sec -2 . 

TRIAL DATA FOR PROGRAM CHECKOUT 

R = 6378145 m 

G = 6.6732E - 11 m 3 kg -1 sec -2 

A 2 = 0.564 kg m -2 

J = 77 

n = 202 

t* = 50000 sec 

x, = 3151529.23 m 

x 2 = 5458608.75 m 


3639072.50 


From external source (CELEST subroutine PRENUT) 

( - 0.8405285753 0.5417623775 0.2289080162E - 02 

- 0.5417605355 -0.8405316908 0.1413662999E - 02 

0.2689913850E- 02 -0.5190827 376E- 04 0.9999963803 

yj = +316648.61 m 

y 2 = - 6290363.38 m 

y 3 = +3647253.32 m 

t** = 208?3333333 

N = (731 +202) = 933 

AT = 5.612148000E - 04 

t* = 0.5787037037 deg 

d = 2.832607926E + 04 
T = 0.7755257840 
h = 2?819922141E + 04 
s = 3?735060861E + 05 
X = 272?8817616 
h - s = - 3?4530686469E + 05 
o = -3?448331496E + 05 
2a = - 6?896662992E + 05 
x = 0.5011240256 


P^(x) = 2.246624133 
P^(x) = 4.256662025 


3U 

dr 

3U 

3X 

3U 

30 

3r 
3yi 
3X 

3y i 

30 

3y l 

3r 

3y 2 

3\ 

ay 2 

30 

3y 2 

3r 

ay 3 

3X 

9y 3 

30 

3y 3 


= 7.070049216E- 12 

= - 5.416307419E - 04 

= 2.46749468IE - 05 

= 4.350677408E - 02 

= 1.585715104E - 07 

= - 3.461599518E- 09 

= - 8.642811299E- 01 

= 7.982281040E- 09 

= 6.876619114E- 08 

= 5.011240258E - 01 

= 0 

= 1.189005543E - 07 

■ - 8.566502457E - 11 


m sec" 


m 2 sec 


m 2 


sec 


8.737156821E - 12 
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T yj = 6.476836379E - 12 m sec -2 

T x = 7.675476994E - 11 m sec" 2 

= 7.675476994E - 14 km sec" 2 

T X2 = - 3.906656638E- 11 m sec" 2 

= - 3.906656638E - 14 km sec" 2 

T x = 6.26836743IE - 12 m sec" 2 

= 6.26836743IE - 15 km sec -2 

T = |f| = 8.635267078E - 11 m sec" 2 

= 8.635267078E- 14 km sec' 2 


2 - COMPUTER ROUTINE FOR THE PERTURBING ACCELERATION 
CAUSED BY THE SOLAR AIR TIDE 


INPUT DATA 


R = any reasonable value for the mean radius of the earth, in meters 
G = gravitational constant, in meters, kilograms, and seconds 
suggested value: G = 6.6732E - 11 m 3 kg " 1 sec -2 
A, I constants associated with the physics of the solar air tide 
? = suggested values: Aj = 6 kg m ~ 2 
A 2 J A 2 = 11.9 kg m ~ 2 

The parameters J, n, and t* specify the time instant for which the perturbing 
acceleration due to the tide is to be found. Normally, this is the time instant associated 
with the current time line of the orbit integration. 

J = number of the calendar year 
n = number of the day within the year 
t* = mean solar time at Greenwich (GMT, UTC), in seconds 

x, 1 inertial, Cartesian components of the satellite position vector, in meters, 
x 2 > = associated with the just specified time instant (for a precise definition 
x 3 j of the reference frame, see Appendices B and C of Reference 2) 
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ALGORITHM 


First, calculate the earth-fixed, Cartesian satellite coordinates (yj, y 2 , and y 3 ) from 
the corresponding inertial coordinates (x,, x 2 , and x 3 ). 


( ABCD )j.n.«-^2^ (M) 


where (ABCD)j n is the transformation matrix that manages the transition from the 
inertial frame of the satellite equations of motion (Basic Inertial System , which is 
associated with either 1950.0 or with the beginning of the day UT of trajectory epoch) to 
the earth-fixed reference frame corresponding to the time instant, J, n, t*. This 
transformation is a computer routine that is already part of TERRA/CELEST. It is external 
to the present algorithm. For details, see Appendices B anc C of Reference 2. 


In the earth-fixed frame, the components of the perturbing acceleration are 


3U 

3r 

3U 3X 

3U 

30 

37 

ay? 

3X dy j 

+ 30 

ay, 

3U 

dr 

3U dX 

3U 

30 

dr 

ay 2 

3X dy 2 

+ 30 

ay 2 

3U 

dr 

3U 3X 

dU 

30 

3r 

ay 3 

3X dy 3 

+ 30 

ay 3 


where 



— = . y2 

d yi y] + y\ 

ax = yi 

ay 2 yf+y 2 


( 2 - 2 ) 

(2-3) 

(2-4) 


(2-5) 

( 2 - 6 ) 

(2-7) 
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ay 3 

= 0 

( 2 - 8 ) 

be 

3Yi 

= ~ y i y3 

r*Vyf + y\ 

(2-9) 

30 

3y 2 

- vi ys 
r 2 Vyf+yf 

( 2 - 10 ) 

be 

3y 3 

VyT+yf 

r 2 

( 2 - 11 ) 

3U 

3r 

4r4 jr3 

= a, -p-P‘(sin0) cosa-a 2 — 4 ~P 2 (sin 0 ) cos 20 

5r5 

+ a 3 -jg- P 2 (sin 0) cos 20 

( 2 - 12 ) 

3U 

3X 

R 4 2R 3 

= a t —j P|(sin 0) sin a - a 2 —j- P^(sin0) sin 20 



2 rJ 

+ a 3 -jj- P 2 (sin 0) sin 20 

(2-13) 

3U 

30 

= a, (15 sin 2 0 - 11) sin 0 cos a - a 2 sin 20 cos 20 

+ a 3 —jj— (7 sin 2 0-4) sin 20 cos 20 

(2-14) 

r 

■ + y 2 +y 3 

(2-15) 

•i 

. 8 jtGR 

= A.- 

1 105 

(2-16) 

a 2 

5jt 2 GR 

= A2_ ir 

(2-17) 

a 3 

_ 5jt 2 GR a 2 

2 3072 " 48 

(2-18) 

Pj(x) 

= | Vi - x 2 (5x 2 - 1) 

(2-19) 
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P^(X) = 

3(1 - x 2 ) 

( 2 - 20 ) 

P^(X) = 

y(l- x 2 )(7x 2 - 1) 

( 2 - 21 ) 

X = 

sin 0 = 

( 2 - 22 ) 

sin 20 = 

2y 3 ^yl+y] 

r 2 

(2-23) 

cos a = 

cos a* cos X - sin a* sin X 

(2-24) 

sin a = 

cos a* sin X + sin a* cos X 

(2-25) 

cos 20 = 

cos 2 0 - sin 2 0 

(2-26) 

sin 20 = 

2 sin 0 cos 0 

(2-27) 

COS0 = 

cos 0* cos X - sin 0* sin X 

(2-28) 

sin 0 = 

sin 0* cos X + cos 0* sin X 

(2-29) 

cosX = 

Vi 

(2 30) 

Vyf+y 2 


sin X = 

V 2 

(2-31) 

Vyf + y 2 

a* = 

t** - 78° 

(2-32) 

0 * = 

t** - 146° 

(2-33) 

*** " 

360 0 

(2-34) 


86400 * 


Finally, transform the perturbing acceleration back into the inertial frame 



(2-35) 
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where (ABCD) T is the transpose of (ABCD). 

While coding the above algorithm for use with TERRA and CELEST, it was found 
helpful to keep in mind the following. A time line is characterized by the time elapsed 
since the epoch (starting instant) of the particular orbital arc. The epoch is specified in 
terms of year, number of the day UT in the year, and time (t EP ), in seconds, elapsed 
from the start of the day UT until epoch time. Usually, t Ep =0. TTie present algorithm is 
periodic in t*, the latter being the Universal time associated with the individual time line, 
t* starts from zero at the beginning of the day UT during which the time line occurs. It 
is thus unimportant how many integer days have elapsed since epoch. For an illustration, 
see Figure 1. 


EPOCH TIME LINE 

I I 



EPOCH integer number of days day lit 

DAY UT 0F 


TIME LINF 

Figure 1. Universal Time in the Algorithm for the 
Perturbing Acceleration Due to the Solar Air Tide 


To actually find the numerical value of t* corresponding to the time line, consider 

that 


t EP = Epoch time (Universal time at epoch) 

At = time difference between time line and epoch 
t* = Universal time of time line 

Now, express t Ep (if non-zero) and At and t* in terms of a common time unit (usually 
mean solar seconds). Find 
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(2-36) 


t' = t £p + At 

Express t' in terms of days and decimal fraction of days (1 
t* = FRACT(t') 


Finally, convert t* to angular time 


t** = 360° FRACT(t') 

which may be substituted for Equation 2-34. 

COMPUTER PROGRAM OUTPUT 

The inertial components (T x j, T x 2 , and T x 3 ) of the 
the presence of the solar air tide are output in terms of n 
km sec -2 . 

TRIAL DATA FOR PROGRAM CHECKOUT 


R = 6378145 
G = 6.6732E - 11 
A, = 6 
A 2 = 11.9 
J = 1977 
n = 202 
t* = 50000 


m 

m 3 kg" 1 sec" 
kg m" 2 
kg m“ 2 


sec 


day = 86400 sec). Now, 

(2-37) 

(2-38) 

perturbing acceleration due to 
i sec” 2 . If required, convert to 
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X. = 3151529.23 


m 


x 2 = 5458608.75 m 

x 3 = 3639072.50 m 


From external source (CELEST subroutine PRENUT) 


( - 0.8405285753 
- 0.5417605355 
0.2689913850E - 02 


0.5417623775 
-0.8405316908 
-0.5190827376E- 04 


0.2289080162E - 02 
0.1413662999E - 02 
0.9999963803 


Vt 


y 2 


y 3 

r 

X 


t** 

a 

0 

Pj(x) 

P^(x) 

P^(X) 

a l 

a 2 


= 316648.61 
= - 6290363.38 
= 3647253.32 
= 7278145.00 
= 272®8817616 
= 30^07439285 
= 208°3333333 
= 403®215095 
= 335®215095 
= 0.3318192840 
= 2.246624133 
= 4.256662025 
= 6.1126614I3E - 04 
= 3.905397704E - 03 
= 8.136245217E- 05 


m 

m 

m 


nr sec d 
m 2 sec ' 
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au 

3r 

3U 

ax 

au 

30 


30 

ay, 

3r 

ay 2 

ax 

ay 2 

it 

ay 2 

3r 

ay 3 

ax 

ay 3 

30 

ay 3 


= - 1.450823033E- 09 

= 8.79906659IE - 03 

= - 6.659222442E - 03 

= 4.350677408E - 02 

= 1.585715104E - 07 

= - 3.461599518E- 09 

= - 8.642811299E- 01 

= 7.982281040E- 09 

= 6.876619116E- 08 

= 5.011240258E - 01 

= 0 

= 1.189005543E - 07 

= 1.355212210E - 09 


m sec 2 
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T y2 = 8.662262286E - 10 m sec " 2 

T y3 = - 1.518827519E - 09 m sec " 2 

T s = - 1.612467289E- 09 m sec " 2 

= - 1.612467289E - 12 km sec " 2 

T X2 = 6.191232115E- 12 msec " 2 

= 6.191232115E- 15 km sec " 2 

T Xj = - 1.514495280E - 09 m sec " 2 

= - 1.514495280E - 12 km sec " 2 

T = |T| = 2.21219010IE - 09 m sec " 1 

= 2.212190101E - 12 km sec " 2 


3 - FIRST COMPUTER ROUTINE FOR THE PERTURBING ACCELERATION 
CAUSED BY THE M 2 COMPONENT OF THE OCEAN TIDE (BASED ON A 
DISCRETE REPRESENTATION OF THE SEA SURFACE) 


INPUT DATA 

fjj = tidal amplitude on area element y (these are output data from the 
NSWC/Schwiderski Ocean Tide Program); if available in meters, use 
directly; otherwise, convert to meters 

8 ^ = tidal phase angle on area element ij (these are output data from the 
NSWC/Schwiderski Ocean Tide Program), available in terms of degrees 

NP = number of grid points in the NSWC/Schwiderski Ocean Tide Program; 
this may be expected to be a Five-digit integer 

R = radius of the earth (semimajor axis of a suitable reference ellipsoid), 
in kilometers 

suggested value; R = 6378.145 km 
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e 2 = square of eccentricity of reference ellipsoid 
suggested value: e 2 = 0.00669342 

in case it is desired to start from the flattening, f, find e 2 from 
e 2 = (2 - Of 

H E = gravitational constant of the earth, in kilometers and seconds 
suggested value: ' 398601 km 3 sec -2 

G = Newton’s constant, in kilometers, kilograms, and seconds 
suggested value: G = 6.6732E-20 km 3 kg -1 sec -2 

p = density of water, in kilograms and kilometers 
suggested value: p = l.E+12 kg km -3 

a = rate of mean mean longitude of moon 

180 . 
suggested value: a = — 1.40519E- 04 deg sec 1 

7r 

The parameters J, n, and t* specify the time instant for which the perturbing 
acceleration due to the tide is to be found. Normally, this is the time instant associated 
with the current time line of the orbit integration 

J = number of the calendar year 
n = number of the day within the year 
t* = mean solar time at Greenwich (GMT, UTC), in seconds 

Xj \ inertial, Cartesian components of satellite position vector, in 

x 2 ( = kilometers, associated with the just specified time instant (for a 
x 3 l precise definition of the reference frame, see Appendices B and C of 
' Reference 2) 

NMAX = limit on the number of terms in the expansion of the tide potential 


SUBROUTINE FOR THE TIDE POTENTIAL COEFFICIENTS 


This is the first part of the algorithm for the first ocean tide computer routine. The 
time-independent constituents of the expansion coefficients of the ocean tide potential 
(“ F nm’ ?F nm> “ H nm > and ^ H nm ^ are calculated from the amplitude and phase tables 
produced by the NSWC/Schwiderski Ocean Tide Computer Program. Once established, the 
latter remain unchanged as long as the current versions of the amplitude and phase tables 
remain valid. Thus, the subroutine under discussion will be exercised quite infrequently, 
namely, just once each time a new edition of the just mentioned amplitude and phase 
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tables becomes available. That is expected to occur at intervals of several years, during 
which the present subroutine will remain dormant. 

For the geometry associated with the index ij, see Figures 2, 3, and 4. 


a (j = 10~ 3 p G ASj 

j fij cos5 ij 

(3-1) 

/J (j = 10~ 3 p G ASj 

j fij si" 5 ,j 

(3-2) 

/ n \ 2 , 

/ 7T \ 


AS.. = — 1 R 

sin { —-j ) 

(3-3) 

,J \180/ 

\ 180 7 

= area of the 

nonpolar surface area element for 


j = 2,3,4 •• 

' .j M AX< 180 


< AS uWar = 2 (lio) 

R 2 

(3-4) 

= area of the 

polar surface area element (neighboring 

the North 

Pole) 

V 7T / 

1 \ 


6. = - - — j 

- - 1 in radians 

(3-5) 

‘J 2 180 \ 

2 ) 

n / 1 \ 

« 180 V" 2) 

| in radians 

(3-6) 

Pii = (l - ^sin 2 

%) 

(3-7) 


To each index ij assign now a counting number v, 1 < v < NP. 


“ F nm 

= (2 . 6 o —L_ 

(2 5,n) (n + m)! p E R 2n 

_ n (n - m)! 1 

N P 

L P 2 v n 

NP , 

+ l a l ,r m 

(3-8) 

fF nm 

= (2 - 8 ®,) ----y- 

m (n + m)! p E R 2n 

(n - m)! 1 " p , 

L p ]" 

+ ' Wm 

(3-9) 

° H nm 

= - - - r- T P 2 ' 

(n + m)!p E R 2n M v 

1+ ‘ a V h nm 


(3-10) 
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ire 3. Position of Point Mass at the Geometrical Center of the 
Surface Area Element Associated With ii 









Figure 4. Position of the Point Mass M 
Earth-Fixed Cartesian Coordinate Fr 


„ H _ (n- m)! 1 ff ^ ^ 

nm (n + m)! ja E R 2n M v v nm 


■! .(' if 8 = k 

10 S * k 

Calculate the f* 

and h* m as follows, noting that 


R 

P = — 

P„ 

and 

hj! =0 for all values of n 






and that 


h’m = f^ m = 0 for any n < m 


(3-15) 


Obtain, separately for each v, the required f£ m and h* m from the following 
recurrence relations: 

To advance in n, evaluate 


f n+l,m = n - rn + 1 [^ 2n + ^ sin f n,m ~ ( n + m ) P„ f n- 1 m] (3*16) 

h n + l,m = [ (2n+ 1)sin0 , h n,m ' (n + m) P, h n-I,m] (3-17) 

To advance in m, use 

f„+i n+i = ( 2 n + 1 )p y |cos 6 v cos \ p P n - cos 0 v sin \ v h* n J (3-18) 

h n+1,n+1 = ( 2 n+ Dp^jcose,, cosX^ h* n + cos 0 p sinX, P n ] (3-19) 

Start these recurrences from 


P„ 



(3-20) 

(3-21) 


(3-22) 


and terminate the procedure at n = NMAX. 

Store the resulting “F nm , ^F nm , “H nm , and ^H nm . They will be programmed as 
constant parameters into the subroutine for the perturbing acceleration and will be expected 
to serve for all subsequent runs of that subroutine, until updated. Note that the “F 
0 F nm , “H nm , and fi H nm are dimensionless quantities. 
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While it would be decidedly impractical to compute F m and h" m from the following 
equation, for program checkout purposes it is possible to find individual, numerical values 
for these quantities from 


/ f nm\ R n 

/ cos mX\ 

( kJ = P "’ (sm0, ' ) 

ysin mX v J 


This relationship may also be used to justify Equation 3-15, considering that 
P™(x) = 0 for n < m 

which is one of the basic identities of Legendre function theory. 


(3-23) 


(3-24) 


SUBROUTINE FOR THE PERTURBING ACCELERATION 

First, calculate the earth-fixed, Cartesian satellite coordinates (y j, y 2 , and y 3 ) from 
the corresponding inertial coordinates (x,, x 2 , and x 3 ). 



(3-51)+ 


where (ABCD)j n t . is the transformation matrix that manages the transition from the 
inertial frame of the satellite equations of motion (Basic Inertial System, which is 
associated with either 1950.0 or with the beginning of the day UT of trajectory epoch) to 
the earth-fixed reference frame corresponding to the time instant, J, n, t*. This 
transformation is a computer routine that is already part of TERRA/CELEST. It is external 
to the present algorithm. For details, see Appendices B and C of Reference 2. 

Find 


r = + Vx]"+~x|~+"xJ" 


(3-52a) 


r = +\/yf+yf+7|' 


(3-52b) 


as convenient 


Numbering sequence advanced 


i emphasize parallelism between certain equations of the algorithms in this and in the 
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Calculate now the time dependent expansion coefficients for the tide potential from 



F nm = “F nm cos (at* +x) + 0 F nm sin (at* + x) 

(3-53) 


H nm = “H nm cos (at* +X) + p H nm sin (at* + x) 

(3-54) 

Obtain x a $ follows 

[ 1 n = m 


5(n,m) = 

if 

(3-55) 


' 0 n ¥= m 


Nj. = 

N 6(J, 1975) +(365 + N)6(J, 1976) + (731 + N) 5(J, 1977) 

+ (1096 + N) 6(J, 1978) + (1461 + N) 6(J, 1979) 

(3-56) 

AT = 

m 

o 

+ 

u> 

o 

00 

(3-57) 

^ - 

27392.5 + Nj. + AT 

(3-58) 


dp 

(3-59) 

T 0 = 

36525 

X = 

270.434358 + 481267.88314137 T 0 - .001133 Tj* + .0000019 T* 

(3-60) 


Note that the day number N, and thus x, must be updated whenever t* runs into the 
next day. 

Remember that F nm and H nm are linear functions of two trigonometric terms. To 
evaluate these trigonometric terms, let t* and t* + t be the values of Universal time for 
which subsequent integration steps (time lines) are to be performed. Update the 
trigonometric time factors as follows: 

cos (at* + , + x) = cos j (at* + x) + oAt } = cos (ot* + x) cos oAt - sin (ot* + x) sin a At (3-61) 
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sin (ot* + ( +x) = sin {(at* + x) + oAt} = sin (at* + x) cos aAt + cos (at* + x) sin aAt 


Now calculate the Cartesian components of the perturbing acceleration 
earth-fixed frame 


d4> 

ay, 

NM AX 

E 

n = 0 

A 

(r.. 

3y, 

av„ 

ay, 

M _ 

NM AX 

n 

I 

aUnm 

av n 

ay 2 

E 

n = 0 

E 

m =0 

l P “ 

~ 12 + H 
ay 2 

ay 2 

M 

NM AX 

V* 

L 

au nm 

av n 

ay 3 

E 

n = 0 

E 

m =0 

V F - 

X- + H m 

ay 3 

ay 3 


r) 


au nn 

ay, 

au nn 

ay 2 

3U nit 

dy 3 


1 / i i_ \ 

1 i 2 ^ n +1 -m - 1 2 + I ,m ♦ 1 ) 

i(-I A V -Iv \ 

• \ 2 m n n+l.m-l 2 V n+l,m + l/ 


av ntn l A \ 

3y, R A n>n V n + i,m-i ' 2 V " + i.« n + ij 

av nm _ i/_i J_ \ 

3y 2 R l 2 AmnU " +1 -" , -i + 2 U " +| .m + | / 


av nr 

»y 3 


~(n - m + 1 )V n , 


(3-62) 

(3-63) 

in the 

(3-64) 

(3-65) 

(3-66) 

(3-67) 

(3-68) 

(3-69) 

(3-70) 

(3-71) 

(3-72) 
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and 


- m + l)(n - 

m + 2 ) 

(3-73) 

(“-»>» » 

Uni 

(3-74) 

(n+ 1 )! nl 

n(n + 1 ) 

Di 

Vnl 

(3-75) 

(n+ 1 )! nl 

n(n + 1 ) 


Obtain the values of the individual eigenfunctions from the following recursive 
relationships: 

(3-76) 
(3-77) 

V no = 0 for all values of n (3-78) 

U nm = V nm =0 for all n < m (3-79) 



To advance in n, evaluate 


u = —-- 

u n + i. m n - m + 1 

|(2n + 1) sin U nm 

- (n + m)p U n _, m [ 

(3-80) 

V = p 

n +1 ,m n- m + 1 

j(2n + 1) sin 0 V nm 

- (n + m)p V fl _ j m | 

(3-81) 

To advance in m, use 




U„ +Itn+1 = (2n + 1 )p | 

V - —V ) 
{ r nn r »»/ 


(3-82) 
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(2n + 1)P^~T V nn + Y - t U nn ) 


(3-83) 


Start from 


Ry _ 3 

' H r 3 


(3-84) 

(3-85) 

(3-86) 


The Cartesian components of the perturbing acceleration, in the earth-fixed reference 
frame, are now 


00 

(3-87) 

3y, 


00 


3y 2 

(3-88) 

00 

(3-89) 

9y 3 



Finally, rotate the perturbing acceleration back into the inertial coordinate system 



(3-90) 


where (ABCD) T is the transpose of (ABCD). 
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COMPUTER PROGRAM OUTPUT 


The inertial components (T x ( , T X2 , and T X l ) of the perturbing acceleration due to 
the presence of the ocean tide are output in terms of km sec' 2 . If required, convert to 


TRIAL DATA FOR PROGRAM CHECKOUT 

Non-zero tidal heights and associated phase values were assumed for nine surface area • 
elements. Several of the “F nm , 0 F nm , “H nm , and tf H nm were subsequently computed. The 
algorithms for the perturbing acceleration, featured in the first and second computer 
routines for the ocean tide, are identical. Trial data for these segments of the two tide 
routines will therefore be listed only once, namely, for the second ocean tide routine, 
below. 


Quantities Associated With the Surface Area Elements 




AS „ 

w 


1.1 I 10 25 

2.1 2 10 25 

3.1 3 10 25 

1.2 4 10 25 

2.2 5 20 30 

3.2 6 20 30 

1.3 7 10 25 

2.3 8 20 30 

3.3 9 20 30 


1.562069681 8.7266463 E 03 

1.562069681 2.61799388E 02 

1.562069681 4.36332313E - 02 

1.544616388 8.7266463 E - 03 

1.544616388 2.61799388E - 02 

1.544616388 4.36332312E - 02 

1.527163095 8.7266463 E - 03 

1.527163095 2.6I799388E 02 

1.527163095 4.36332313E - 02 


6356.800824 108.1411251 

6356.800824 108.1411251 

6356.800824 108.1411251 

6356.813825 432 4766612 

6356.813825 432.4766612 

6356.813825 432.4766612 

6356.839812 648.550316 

6356.839812 648.550316 

6356.839812 648.S50316 


6.5403462E - 08 
6.5403462E - 08 
6.5403462E - 08 
2.6156072E - 07 
4.9987043E - 07 
4.9987043E - 07 
3.9224149E -07 
7.496153 E - 07 
7.496153 E- 07 


3.0498135E - 08 
3.0498135E- 08 
3.0498135E - 08 
1.2196777E - 07 
2.8860033E - 07 
2.8860033E - 07 
1.8290S21E - 07 
4.327906 E - 07 
4.327906 E - 07 





Values for the f* and h" 


V 

**0,0 
km 1 

h o,o 
km 1 

km 1 

,0 

km 1 

l 

1.5731184E - 04 

0 

1.5783403E - 04 

0 

2 

1.573II84E- 04 

0 

1.5783403E - 04 

0 

3 

1.5731184E-04 

0 

1.5783403E - 04 

0 

4 

l .5731 15 IE - 04 

0 

1.5778531E-04 

0 

5 

1.573115 IE - 04 

0 

1.577853IE - 04 

0 

6 

1.573115 IE - 04 

0 

1.577853 IE-04 

0 

7 

1.5731087E-04 

0 

1.5768788E - 04 

0 

8 

1.5731087E - 04 

0 

1.5768788E - 04 

0 

9 

I.5731087E - 04 

0 

I.5768788E-04 

0 


V 

*1.1 
km " 1 

h\ , 
km' 1 

* 2.0 
km" 1 

h 2,0 
km 1 

1 

1.3773442E - 06 

1.2019901E-08 

1.5835193E - 04 

0 

2 

1.3769246E - 06 

3.605604 1 E - 08 

1.5835193E - 04 

0 

3 

1.3760857E-06 

6.0081198E-08 

1.5835193E - 04 

0 

4 

4.1315963E - 06 

3.6055895E - 08 

I.5820627E-04 

0 

5 

4.1303378E - 06 

I.0815670E -07 

1.5820627E - 04 

0 

6 

4.127821 IE - 06 

1.8022456E - 07 

1.5820627E - 04 

0 

7 

6.8845393E - 06 

6.0080464E - 08 

I.57915I3E - 04 

0 

8 

6.8824422E - 06 

1.8022309E - 07 

1.5791513E - 04 

0 

9 

6.8782468E - 06 

3.0031082E-07 

1.5791513E - 04 

0 


„ 

*2,1 
km" 1 

h 2 , 1 
km 1 

f 2,2 
km 1 

h 2.2 
km 1 

1 

4.1457488E - 06 

3.6179402E - 08 

3.6175267E - 08 

6.3I44163E- 10 

2 

4.1444859E - 06 

1.0852718E-07 

3.6131193E-08 

1.893 5 5 56E - 09 

3 

4.1419607E - 06 

1.8084191E -07 

3.6043099E - 08 

3.I533625F. - 09 

4 

1.2432120E - 05 

1.0849347E - 07 

3.2550933E - 07 

5.68I7864E -09 

5 

1.2428333E - 05 

3.2544735E - 07 

3.2511274E - 07 

l .7038437E - 08 

6 

1.2420760E - 05 

5.4230210E-07 

3.2432006E - 07 

2.8374329E - 08 

7 

2.0703116E-05 

1.8067335E - 07 

9.0381431 E-07 

1.5776137E - 08 

8 

2.0696809E - 05 

5.4I96503E-07 

9.0271315E - 07 

4.730919IE - 08 

9 

2.0684198E-05 

9.0309161 E-07 

9.00512I7E-07 

7.8784606E - 08 
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Values for the fJJ m and h* m 

(Continued) 



f 3,0 

h 3,0 

f 3',. 

h 3 

V 

km' 1 

km - 1 

km” 1 

km 


1 

1.5886548E-04 

0 

8.3188628E - 06 

7.2597615E-08 

2 

I.5886548E-04 

0 

8.3163287E - 06 

2.1777073E - 07 

3 

I.5886548E-04 

0 

8.3112614E - 06 

3.628775 IE-07 

4 

1.5857391E - 04 

0 

2.4934851E - 05 

2.1760315E - 07 

5 

1.585739IE - 04 

0 

2.4927255E - 05 

6.5274316E-07 

6 

1.5857391E - 04 

0 

2.4912067E - 05 

1.0876843E - 06 

7 

1.5799154E - 04 

0 

4.1485684E - 05 

3.6204008E - 07 

8 

1.5799154E-04 

0 

4.1473047E - 05 

1.0860100E - 06 

9 

1.5799154E - 04 

0 

4.1447777E - 05 

1.8096490E - 06 


V 

f 3,2 

km 1 

h 3.2 

km 1 

f 3,3 
km 1 

h $,3 
km -1 

1 

1.8147675E-07 

3.1676885E - 09 

1.5834220E - 09 

4.1463365E- 11 

2 

1.8125565E - 07 

9.4992060E - 09 

1.579082 E - 09 

1.2427645E - 10 

3 

1.8081372E - 07 

1.5819151E - 08 

1.5704138E - 09 

2.0674889E - 10 

4 

1.6324485E - 06 

2.8494495E - 08 

4.2739029E - 08 

1.1191609E - 09 

5 

1.6304596E-06 

8.5448767E - 08 

4.2621885E - 08 

3.354415IE - 09 

6 

1.6264843E - 06 

I.4229893E - 07 

4.2387916E - 08 

5.58047 E - 09 

7 

4.5299018E - 06 

7.9069730E - 08 

l 9774213E - 07 

5.1780599E - 09 

8 

4.5243828E - 06 

2.371 I286E-07 

1.9720013E- 07 

1.5519987E - 08 

9 

4.5133516E - 06 

3.9486710E-07 

1.9611762E-07 

2.5819375E - 08 


V 

f ;.3 
km" 1 

h S,3 
km 1 

1 

1.1120747E - 08 

2.912070IE - 10 

2 

1.1090266E - 08 

8.7282284E - 10 

3 

1.1029387E - 08 

1.4 5 20463 E - 09 

4 

3.0007426E - 07 

7.857721 IE-09 

5 

2.9925178E-07 

2.3551626E - 08 

6 

2.9760907E - 07 

3.9180977E - 08 

7 

1.3875122E-06 

3.6333286E - 08 

8 

1.383709IE - 06 

1.0890027E - 07 

9 

1.3761134E-06 

1.8116877E - 07 
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Time-Independent Coefficients of the Tide Potential 


n 

m 

“F nm 




a » nm 




0 

0 

8.4018456E- 

12 

4.6140106E- 

12 

0 


0 


1 

0 

8.368164IE - 

12 

4.5954868E - 

12 

0 


0 


1 

1 

2.9297877E - 

13 

1.6202500E - 

13 

4.3123569E- 

15 

2.4544968E - 

15 

2 

0 

8.3290386E- 

12 

4.5739464E - 

12 

0 


0 


2 

1 

2.9177584E - 

13 

1.6I35948E - 

13 

4.2946458E - 

15 

2.4444118E — 

15 

2 

2 

2.7840515E - 

15 

1.5423222E- 

15 

8.2132887E - 

17 

4.683368 E- 

17 

3 

0 

8.2845357E - 

12 

4.5494263E - 

12 

0 


0 


3 

1 

2.904663 2E- 

13 

1.6063487E - 

13 

4.2753632E - 

15 

2.4334303E - 

15 

3 

2 

2.7724443E - 

15 

1.5358913E — 

15 

8.1790436E - 

17 

4.6638388E - 

17 

3 

3 

1.581308IE - 

17 

1.0259185E - 

17 

8.2068535E - 

19 

4.6816230E - 

19 

4 

3 

1.8435105E- 

17 

1.0215973E - 

17 

8.1722863E - 

19 

4.66I9036E- 

19 


4 - SECOND COMPUTER ROUTINE FOR THE PERTURBING ACCELERATION 
CAUSED BY THE M 2 COMPONENT OF THE OCEAN TIDE (BASED ON A 
FUNCTIONALIZATION OF THE SEA SURFACE) 


INPUT DATA 


C„ I 

Sjj I NOAA/Goad expansion coefficients of the tidal surface, in meters; 
C,j j these are specified by a computer card deck 

S 9 ) 

R = radius uf the earth (semimajor axis of a suitable reference ellipsoid) 
in kilometers 

suggested value: R = 6378.145 km 


G = Newton’s constant, in kilometers, kilograms, and seconds 
suggested value: C. = 6.6732E 20 km 3 kg ' 1 sec ' 2 


Hy = gravitational constant of the earth, in kilometers and seconds 
suggested value: n v = 398601 km 3 sec 2 

p = density of water, in kilograms and kilometers 
suggested value: p = l.E+12 kg km 3 
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a = rate of mean mean longitude of moon 
180 

suggested value, o = — 1.40519E- 04 deg sec 1 

7T 

NMAX = limit on the number of terms in the expansion of the tide potential; 

NMAX is equal to the largest value of the index i for the expansion 
coefficients of the tidal surface 

The parameters J, n, and t* specify the time instant for which the perturbing 
acceleration due to the tide is to be found. Normally, this is the time instant associated 
with the current time line of the orbit integration. 

J = number of the calendar year 
n = number of the day within the year 
t* = mean solar time at Greenwich (GMT, UTC), in seconds 

x,| inertial, Cartesian components of satellite position vector, in 

x 2 l = kilometers, associated with the just specified time instant (for a 
x 3 I precise definition of the reference frame, see Appendices B and C of 
' Reference 2) 


SUBROUTINE FOR THE TIDE POTENTIAL COEFFICIENTS 

This is the first part of the algorithm for the second ocean tide computer routine. 
The time-independent constituents of the expansion coefficients of the ocean tide potential 
(Fjj, F", H'., H'j) are calculated from the expansion coefficients for the NOAA/Goad tidal 
ocean surface. Once established, the latter remain unchanged as long as the current version 
of the NSWC/Schwiderski Ocean Tide Model (on which the NOAA/Goad tidal surface is 
based) remains valid. Thus, the subroutine under discussion will be exercised quite 
infrequently, namely, just once each time a new edition of the NSWC/Schwiderski Ocean 
Tide Model becomes available. That is expected to occur at intervals of several years, 
during which the present subroutine will remain dormant. 
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H " = 
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2 
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SUBROUTINE FOR THE PERTURBING ACCELERATION 

Proceed exactly as in Equations 3-51 through 3-90 from the first ocean tide computer 
routine, but delete Equations 3-53 and 3-54 and replace them with 


= F' m cos (at* + x) + F” m 

sin (at* + x) 

(4-53) 

= H' nm cos (at* +x) + H' n ' ni 

, sin (at* +x) 

(4-54) 


COMPUTER PROGRAM OUTPUT 

The inertial components (T x f , T X2> and T Xj ) of the perturbing acceleration due to 
the presence of the ocean tide are output in terms of km sec -2 . If required, convert to 
m sec -2 . 


TRIAL DATA FOR PROGRAM CHECKOUT 

Non-zero values were assigned to selected NOAA/Goad expansion coefficients, and the 
corresponding expansion coefficients of the tide potential were computed. The algorithm for 
the perturbing acceleration was then exercised with the following results. 
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50000 

sec 


6378.145 

km 


6.6732E - 20 

km 3 

kg - 

398601 

km 3 

sec' 
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p = 10' 


kg km 3 

180 

a = -1.40519E - 04 deg sec~ 

NMAX = 4 

C 20 = + 0.2906060089E - 01 m 
C 40 = - 0.107121752E + 00 m 
C 43 = + 0.435761219E - 04 m 

C 2 ' 0 = - 0.4424413130E-01 m 
C 40 = + 0.873468034E - 01 m 
C 43 = - 0.160563906E-02 m 

s 20 - 0 

S 40 = 0 

S 43 = - 0.363303008E - 02 m 
S 2 ’ 0 = 0 

S 4 'o = 0 

s 4 ' 3 = - 0.264356490E-02 m 

F 20 = + 4.9742658E - 10 
F 40 = - 1.0186607E - 09 
F 43 = +4.1438160E- 13 
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r 

Fjto “ - 7.5732111 IE- 10 
F" 0 = + 8.30613340E- 10 
F" 3 = - 1.526862 IE- 11 

L 

H' 20 - 0 

h; 0 - o 

H; 3 = - 3.4547839E - 11 
H''o = 0 

h ;' 0 = 0 

H” 3 = - 2.5138645E - 11 

N £ = 933 
AT = 0.000561 
d Q = 28325.50056 
T 0 = 0.7755099401 
X = 373498.4609 

o = 0.0080511456 deg/sec 

at* = 402.557282 deg 

cos (at* + x ) = -0.754501 
sin (at* + x) = - 0.656298 

F 2q = + 1.2171968E - 10 

h 20 = 0 

F 40 = + 2.2345060E - 10 
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9.7081216E - 12 


From external 


(ABCD) = 


H 43 = + 4.2564846E - 11 


= 3151.52923 km 

x 2 = 5458.60875 km 


x 3 = 3639.07250 


km 


source 


-0.8405285753 
-0.5417605355 
0.2689913850E- 02 


0.5417623775 0.2289080162E - 02 

- 0.8405316908 0.1413662999E- 02 

- 0.5190827376E - 04 0.9999963803 


y, = + 316.64861 km 

y 2 ■ - 6290.36338 km 

y 3 = +3647.25332 km 

dU 2(J 

r —— = - 0.0000964047 km/sec 2 

ay, 

9V 20 

— = - 4.7E - 13 km/sec 2 

du 20 

r —— = +0.0019151219 km/sec 2 

ay 2 

r—— = - 2.4E - 14 km/sec 2 

dy 2 
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Values for the Eigenfunctions of the Tide 
Potential, in km 2 sec -2 


n 


0-1 0 
0 0 54.76683963 


-1.044042677 

24.05119113 

2.088085354 


2 

2 

2 

2 


2 


-0.4584976993 

-5.186455141 

2.750986196 

-94.01442066 


3 

3 

3 

3 

3 


2 

3 


-0.0512402699 

-16.10992273 

0.6148832385 

-206.435027 

-53.85812191 


4 

4 

4 

4 

4 

4 


-1 0.109342534 

0 -9.393545781 

1 -2.18685068 

2 -136.7982658 

3 -165.56486 

4 1863.67849 


5 

5 

5 

5 

5 

5 


2 

3 

4 


0.0917032956 

2.472201758 

-2.751098868 

136.846714 

-182.4236505 

7366.011828 


6 

6 

6 

6 

6 

6 


2 

3 

4 


0.0153004451 

8.002095405 

-0.6426186943 

349.1179421 

45.32033435 

11350.89177 


V n, 


0 

0 

-20.74036524 

0 

-41.48073047 

-9.108257692 

0 

-54.64954615 

-9.489169689 

-1.017910413 

0 

-12.21492495 

-20.83613333 

354.226451 

2.172137348 

0 

43.44274695 

-13.80747707 

1088.924921 

380.085929 

1.821726148 

0 

54.65178443 

13.8123671 

1199.805712 

1502.253452 

0.303950046 

0 

12.76590193 

35.23756643 

-298.0731712 

2314.94556 
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km/sec 2 

km/sec 2 

km/sec 2 

km/sec 2 

km/sec 2 

km/sec 2 

km/sec 2 

km/sec 2 

km/sec 2 


38 




= +0.0572027292 km/sec 2 

ay 3 

0V 43 

r—~ =- 0.3762240313 km/sec 2 

ay 3 

T yi - - 9.632495E- 12 km/sec 2 

T y2 = + 2.443056E - 11 km/snc 2 

T yj = - 1.496932 IE - 11 km/sec 2 

T x = - 5.179392E - 12 km/sec 2 

T Xj = - 2.5752406E- 11 km/sec 2 

T Xj = -1.495676E- 11 km/sec 2 


5 - COMPUTER ROUTINE FOR THE PERTURBING ACCELERATION 
CAUSED BY THE TIDES OF THE SOLID EARTH 


INPUT DATA 

The parameters J, n, and t* specify the time instant t for which the perturbing 
acceleration due to the tide is to be found. Normally, this is the time instant associated 
with the current time line of the orbit integration. 

J = number of the calendar year 
n = number of the day within the year 
t* = mean solar time at Greenwich (GMT, UTC), in seconds 

At = tidal lag parameter, in seconds 

£3 = earth’s sidereal rate of rotation 

suggested value: £> = 4.178074622E - 03 deg sec' 1 
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the inertial, Cartesian coordinates of the moon, in kilometers, for the 
time instant (t - At), from the TERRA/CELEST moon ephemeris; the 
reference frame is referred to mean equator and mean equinox at time 
instant t, 

the inertial, Cartesian coordinates of the sun, in kilometers, for the time 
instant (t - At), from the TERRA/CELEST sun ephemeris; the reference 
frame is referred to mean equator and mean equinox at time instant tj 

t, = epoch time for the moon ephemeris and sun ephemeris, in Julian Date 
t, = JD2433282.423 

x, ) inertial, Cartesian components of satellite position vector, in kilometers, 
x 2 > = associated with the time instant t(J, n, t*) (for a precise definition of 
x 3 \ the reference frame, see Appendices B and C of Reference 2) 

• 1 I _ Cartesian satellite velocity components associated with time instant t, in 
. 2 ( km sec -1 


= gravitational constant of the earth, in kilometers and seconds 
suggested value: /jl £ = 398601 km 3 sec -2 

R = radius of the earth (semimajor axis of a suitable reference ellipsoid) in 
kilometers 

suggested value: R = 6378.145 km 

e 2 = square of eccentricity of reference ellipsoid 
suggested value: e 2 = .006693421623 

in case it is desired to start from the flattening, f, find e 2 from 
e 2 = (2- f)f 

k (j = expansion coefficients for the Love numbers 
suggested values: k 2Q = 0.3 
k 21 = 0.01 
k 30 =0.1 
k 3 , = 0.01 
k 22 = 0.1 
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m' = mass of moon or sun, in kilograms 

suggested values: m ' L = mass of moon = 7.369328IE + 22 kg 
m^ = mass of sun = 1.99E + 30 kg 

M = mass of earth, in kilograms 

suggested value: M = 5.9731613E + 24 kg 

a' = semimajor axis of moon or sun orbit relative to earth, in kilometers; 
select values that are compatible with the moon ephemeris or sun 
ephemeris (optionally, values may be obtained from a current edition of 
the AMERICAN EPHEMERIS ) 
suggested values: a ' L = 384400 km 

a' s = 149500000 km 


SUBROUTINE FOR THE POSITION OF THE FICTITIOUS MOON 



X x = cos £ 0 cos 0 cos z - sin £ 0 sin z 
Y x = - sin £ 0 cos 0 cos z - cos £ 0 sin z 
Z x = - sin 0 cos z 

X = cos £ 0 cos d sin z + sin £ 0 cos z 
Y y = - sin £ 0 cos 0 sin z + cos ,{- 0 cos z 
Z * - sin 0 sin z 
X z = cos sin 6 
Y z = - sin £ 0 sin 0 
Z z = cos0 


(5-1) 

(5-2) 

(5-3) 

(54) 

(5-5) 

(5-6) 

(5-7) 

(5-8) 

(5-9) 

(5-10) 
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* 0 = J(2304.250 + 1.396T 0 )T + 0.302T 2 +0.018T 3 j (5-11) 

2 - t " + ! a79IT2 i 3 -^ (5 ' 12) 

6 = j (2004.682 - 0.853T 0 )T- 0.426T 2 - 0.042T 3 j (5-13) 


£ 0 , z, and 0 result in degrees. 


T 0 = (t, - 2415020.313X0.273790926497E- 04) 
T = (0.273790926497E - 04)(t 2 - t,) 

(t 2 - t,) = j (J - 50)365 + A 0 + (n - 1) + .077 j 




(5-14) 

(5-15) 

(5-16) 

(5-17) 


The Cartesian coordinates of the fictitious moon (x£, y£, and z£) referred to the 
coordinate system in which the satellite orbit integration is performed may now be 
computed: 


/x*\ /cose 

U/Vo 


coscoAt -sin wAt 

coAt cos to At 0 J( y L I 


Also, compute 


r* = +n/(x*) 2 + (y* ) 2 + (z*) 2 
X* = *1 


(5-19) 

(5-20) 


* y* 


(5-22) 





SUBROUTINE FOR THE POSITION OF THE FICTITIOUS SUN 



(5-23) 


Evaluate the matrix elements X x , Y x , etc., according to Equations 5-2 through 5-17. 
Then, calculate the Cartesian coordinates of the fictitious sun (x*, y*, z*) referred to the 
coordinate system in which the satellite orbit integration is performed: 



cos wAt -sinoJAt °\/ x ’s\ 

sin ouAt coscuAt 0 )[ y' s 1 

0 0 \/\z'J 


(5-24) 


Further, compute 


r* = +V(x*)2+(y*)2+( z *)2 



SUBROUTINE FOR SATELLITE POSITION 

r = +Vxf+ r x|TxJ 


(5-25) 

(5-26) 

(5-27) 

(5-28) 


(5-29) 

(5-30) 
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(5-31) 



^3 

r 

(5-32) 

(X,) 2 +(X 2 ) 2 + (Xj) 2 

(5-33) 

/ 2 V 2 

V.-J 

(5-34) 


SUBROUTINE FOR THE PERTURBING ACCELERATION 

T x j, T X2 , and T Xj are the Cartesian components of the perturbing acceleration, in 
km sec- 2 . These are referred to the (inertial) reference frame in which the satellite orbit 

integration is performed, and are thus readily suited for insertion into the TERRA/CELEST 

equations of motion. 

3C, V, C, , C 3 C 4 1 

T x, = 7 )L c o + “T" 1 + tt( 5V 2- 2p o) + tt(7V 3 - F> + ^(9V4-H)Jx, 

rc, , C 2 C 3 C 4 1 I 

‘ lT A3 + 7 P " + 7 rP > 2 + r TP|3 J 1 (5 ' 35) 

-Mf 3C,V, c 2 , C 3 c 4 "I 

T x 2 = 7 j [ c ° + ~i~ A + r T(5V 2 -2P; )+ jj-(7V 3 - F) + ^-(9V 4 -H)Jx 2 

rc, , C 2 C 3 C 4 1( 

" L7 A4 + 7 Pz ' + 7 P22 + 7 p23 J( (5-36) 

T x3 = ~?\[ C o + ^T" 1 + r T-(5V 2 -2P^) + ^-(7V 3 -F)+ ^-(9V 4 -H) x 3 ] 

f- 4C, , C 2 C 3 C 4 "1 ) 

■ |_T~ A o + 7 p 3 . + 7^32 + r — P 33 J j (5-37) 
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2XP' + hiPj + ^P'j 

(1 - 5t> 2 )S' + 3(X 2 - n 2 )S^ + 2X(3/jS' s + i>S^ ) + #u>S 7 
v (' - j * 2 ) r; + 3va 2 - M 2 jt; + 6Xmvt; + a - iv 2 h2xt; + </r' 7 ) 

-2nP\ +XP ' 2 + Kp' 

(1 - 5v 2 )S' - 6 XpiS« + 3(X 2 - ju 2 )S ’ 5 + o(XS 7 - 2nS' b ) 

j*^r 2 - 6 x^i^t; +3wx 2 - // 2 )t; + (i - 7, 2 kxt; - 2pt;> 

- 61 ^ +XP 2 +/ip' 

- lo^xs; + uS' 2 ) + (3 - 15i^ 2 »Sj + (X 2 - » 2 )S' 6 +\fiS' 7 

< 1 - 7t/2 )(XT' t + mT'j ) - 20K3 - lv 2 )T'j 
+ X(X 2 - 3*i 2 )T; +M(3X 2 - *i 2 )T; 

I4i/[(X 2 - H 2 )T^ +X/iT' 7 | 

2 (xs; + j/s; + 3vs'j > 

2 [mxt; + ^t; ) + 6< i sv 2 >t; + (x 2 - n 2 >r + x/jt;] 






and 



V o = 1 

V, = -(AjX + A 'qH 4A| ) i') 
v 2 = 1 P^P n 

v 3 = i^ s^s n 
v 4 = 2 T^T n 

A' = 7 ( 1 - 3t>* 2 ) 

0 4 

A' = 7 (X * 2 n**) 

1 4 

a; = 3X*m* 

a; = 

A^ = 3ij*v* 


(5-54) 

(5-55) 

(5-56) 

(5-57) 

(5-58) 

(5-59) 

(5-60) 

(5-61) 

(5-62) 

(5-63) 

(5-64) 

(5-65) 

(5-66) 

(5-67) 

(5-68) 
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and 


b; 

B' 

K 

b; 

K 

B 1 

B 7 

P 0 

P 0 

P' 

P, 

p ; 

p 2 

P 3 

P 3 

^4 

P 4 


-X*(l - 5k* 2 ) 



(5-69) 

£M*(1 ~ 5k* 2 ) 



(5-70) 

1 

-k*(3 - 5k* 2 ) 



(5-71) 

^*<X* 2 - 3n* 2 ) 



(5-72) 

£M*(3X* 2 - m* 2 ) 



(5-73) 

•5 

— k*(X* 2 - m* 2 ) 

4 



(5-74) 

15X*m*k* 



(5-75) 

( k - t ? k -)( 1 ' 

” £ . 
42 , 


(5-76) 

1 - 3k 2 



(5-77) 

( k >»‘ f k *>)( 1 - 

14 , 

)a; + ^ 31 a'^* B ; 

(5-78) 

X 2 - n 2 



(5-79) 

( k »»- f *»)(■- 

5 > 
14 e > 

)a' 2 + ~ k 3 , a'/3*B' 7 

(5-80) 

Xm 



(5-81) 

( k - k 7 lk »)('- 

£•* 

)a^ - -k 31 a'0* B ' ( 

(5-82) 

Xk 



(5-83) 

( k ’» + ? k »)('- 

15 ,> 

14 f j 

K- “ k ,,«'0*B; 

(5-84) 

(JLV 



(5-85) 
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and 


S' =- jk 2| A' 3 +k 3O a'0*B' | (5-86) 

S, = X(1 - 5f 2 ) (5-87) 

S 2 = - J k 2. A 4 +k 3 0 a '< 3 * B 2 (5-88) 

5 2 = ti( 1 - 5j> 2 ) ( 5 . 89) 

5 3 = ~ k 21 A 0 +k 30 a '< 3 * B 3 (5-90) 

5 3 = v(3 - 5v 2 ) ( 5 - 91 ) 

5 4 = k 30 a ^* B 4 (5-92) 

S 4 = X(X 2 - 3/i 2 ) (5-93) 

s 5 = k 3o a '^* B 5 (5-94) 

S s = n(3\ 2 - n 2 ) (5-95) 

S 6 = k 2. A '. +k 3o a '^* B 6 (5-96) 

S 6 = i>(X 2 - m 2 ) (5-97) 

S ? = k 2i A 2 +k 3o a '0* B '7 (5-98) 

S, = \nv (5-99) 
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(5-100) 

(5-101) 

(50102) 

(5-103) 

(5-104) 

(5-105) 

(5-106) 

(5-107) 

(5-108) 

(5-109) 

(5-110) 

(5-111) 

(5-112) 

(5-113) 


Note that the components T x of the perturbing acceleration depend, essentially, on 
the parameters 


X*, n*. v*, r*. a', m', x,, x 2 , x 3 
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T x = T X .(A*, n*, v *, r*, a', m', 


: i* X 2’ x 3> 


To obtain the perturbing acceleration due to the moon tide, substitute for these parameters 
those specifying lunar position, orbital semimajor axis, and mass. 


(T x .) l = T X .(A£, ii*, v*, r*, a' L , m' L , x,, x 2 , x 3 ) 

Otherwise, substitute the solar parameters to produce the perturbing acceleration caused by 
the tidal action of the sun. 


(T x .) s = T X .(A£, Ms* V S' r s ’ a S’ m S’ X 1 ’ x 2’ x 3) 

COMPUTER PROGRAM OUTPUT 

The inertial components (T X[ ) L , (T X2 ) L , (T X;J ) L , (T X] ) S , (T X2 ) S , and (T x 3 ) s of the 
perturbing acceleration due to lunar and solar tides of the solid earth are input in terms of 
km sec -2 . 


TRIAL DATA FOR PROGRAM CHECKOUT 

The following trial case is restricted to the lunar tide. While executing this case, the 
entire algorithm was exercised and it was, thus, unnecessary to also consider the solar tide. 


Input Data 

t ••• J = 1977 
n = 88 

t* = 57000 sec 

At = 100 sec 

£ = 4.178074622E- 03 deg sec’ 1 


50 




Epoch for the moon ephemeris (reference frame 2,) 
t, = 1950.0 = JD2433282.423 
(1950.0 means beginning of Besselian year 1950) 

Epoch for the satellite orbit (reference frame S 2 ) 


V J = 1977 
n = 88 

t* = 0 sec 

t 2 = JD2443231.5 

Moon coordinates, in frame E,, at time (t - At) 
x" = - 184258.5823 km 

y'i = 329791.1351 km 

z'l = 103840.2349 km 

Satellite position and velocity, in frame Z 2 , at time t 
x, = - 4009.582237 km 

x 2 = + 103.9008135 km 

x 3 = - 5269.570696 km 

Xj = - 5.978210289 km sec' 

x 2 = + 1.710442216 km sec” 

x 3 = + 4.584971066 km sec” 

H E = 398601 km 3 sec 

R = 6378.145 km 


e 2 = 6.693421623E- 03 


k 20 = 0.3 

k 2 , = 0.01 

k 30 = 0.1 

k 3 , = 0.01 

k 22 = 0.1 

m' L = 7.3693281 E +22 kg 

M = 5.9731613E + 24 kg 

a' L = 384400 km 
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Results 


9 = 0.1516444800 deg 

{ = 0.1744119454 deg 

z = 0.1744282487 deg 

T 0 = 0.5000000017 centuries 

T = 0.2723967010 centuries 


\ = +0.9999779632 
Y x = - 0.0060883617 
Z x = - 0.0026466801 
X y = +0.0060883617 
Y y = +0.9999814657 
Z y = - 0.0000080574 
\ = +0.0026466801 
Y z = - 0.0000080567 
Z z = +0.9999964975 


= - 186537.2414 1cm 

y' L = +328662.3531 km 

z' L = + 103349.5407 km 

r* = 391785.9267 km 
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X* = - 188928.9046 km 

y* = +327293.3757 km 

z* = + 103349.5407 km 

X* = - 0.4822248369 

H* = +0.8353882910 

v * = +0.2637908451 

X = - 0.6054593779 

H = +0.0156893457 

v = - 0.7957215507 


r = 6622.380268 km 

a = 6567.451187 km 

K = 0.323073882E - 05 km 2 sec’ 2 

n = 0.118624368E-02 sec" ‘ 


a = 0.1659246878E - 01 
a = 0.9711750904 
0* = 0.9811480551 

A' 0 = +0.1978107925 
A', = - 0.3489996026 
A' 2 = - 0.1208534947E + 01 
Aj = - 0.3816194918 
A; = +0.6611033498 


53 






B', = - 0.1179169837 
B' 2 = +0.2042749770 
B' 3 = +0.1748980753 
B; = +0.5609118737 
B^ = - 0.0001311647178 
B; = - 0.4603145005 
B; = - 0.1594002275E + 01 

Po = +0.06443748370 
P 0 = - 0.5995183587 
P', = - 0.09451271981 
Pj = +0.3663349027 
P' 2 = - 0.3272838250 
P 2 = - 0.009499261487 
P' 3 = - 0.1190554808 
P 3 = +0.4817770751 
p; = +0.2062472668 
P 4 = - 0.01248435049 

S\ = +0.0005712740433 
S, = + 0.1311342628E +01 
S 2 = - 0.0009896538092 
S 2 = - 0.03398098796 
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+ 0.001471593023 


+ 0.1319815046 
+ 0.0009131459348 

- 0.2215028279 

- 0.2135318121E - 06 
+ 0.0172503888 

- 0.004239372772 

- 0.2915005769 

- 0.01468033233 
+ 0.007558767081 

+ 0.02363141132 

- 0.2300019019 

- 0.04093817407 
+ 0.005960068473 

- 0.002454126571 

- 0.1963411384E +01 
+ 0.9131459348E - 04 
+ 0.1762545737 

- 0.213531812IE - 07 

- 0.01372650615 
+ 0.02595922645 

- 0.1257338I35E + 01 
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T’ = +0.08989303178 
T ? = +0.03260345555 

V Q = + 1.0000000000 
V, = +0.1742073243 
V 2 = - 0.1494101174 
V 3 = +0.18995 34786E- 02 
V 4 = - 0.03055341082 

C Q = - 0.0001622651725 
C, = + 0.1353296915E + 01 
C 2 = + 0.8631523953E +06 
C 3 = + 0.5505311134E+ 10 
C 4 = +0.3511367268E+ 14 

F = - 0.007748690189 
H = +0.1048876747 

p,, = +0.2040473678 
p I2 = - 0.004135328808 
p 13 = +0.1119464221 
p 21 = +0.03700735158 
p 22 = - 0.004983233113 


km 3 sec 2 
km 4 sec~ 2 
km 5 sec - 2 
km 6 sec - 2 
km 7 sec 2 
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p 23 = +0.1740434693 
p 31 = +0.3829649087 
p 32 = - 0.01385120787 
p 33 = +0.2036553308 

T x = - 0.1929747594E-09 km sec" 2 

T X2 = + 0.9604765343E - 10 km sec" 2 

T Xj = - 0.1824727285E-09 km sec’ 2 

T = 0.2824193799E- 09 km sec' 2 


DEVELOPMENT STATUS OF THE COMPUTER ROUTINES 


By early summer 1978, all five tide routines had been coded and made part, as 
operational subroutines, of the TERRA and CELEST computer program systems. Program 
checkout was accomplished (separately for each tide routine) by first selecting typical input 
data that would exercise the algorithm to a sufficient degree. Subsequently, the programmer 
would run the algorithm under test on the computer, while the author would perform 
identical calculations on a programmable electronic calculator. Agreement to ten significant 
digits of the computer results with those of the manual computation was taken as evidence 
for successful computer implementation of the algorithm. 

As indicated by the last paragraph of the INTRODUCTION, the two ocean tide 
algorithms may be expected to produce nearly identical perturbing accelerations because the 
tidal ocean surface of the second algorithm is a functionalization of the discrete 
representation of the ocean surface on which the first algorithm is based. Like the routines 
for the two air tides and that for the solid earth tide, each of the two ocean tides had 
been checked out individually. In the process, to facilitate the manual calculations, a simple 
ocean surface had to be selected for each of the two algorithms. Because of the different 
mathematical structure of the two algorithms, simplicity meant that each algorithm required 
its own individual surface. Thus, it was impossible to directly compare the resulting 
perturbing accelerations. This in turn meant that, during this phase of the checkout, the 
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mutual consistency of the two ocean tide routines could not be verified. Because the two 
routines are so complex, this remained a source of concern despite the caution practiced 
while deriving them. The question was finally settled by running each algorithm, using the 
complete input data base. The resulting perturbing accelerations agreed with each other to 
20 percent. Considering the differences in the physical definitions of the tidal mass layers 
(one being a discrete point set, and the other being a continuous surface), this was 
interpreted as satisfactory evidence for the validity of the algorithms. 

For operational use. that ocean tide routine was selected which employs the discrete 
tide surface representation. Realizing that the second algorithm contains a more accurate 
physical model (by virtue of its continuous tidal mass layer as compared with the point 
mass approximation inherent in the first ocean tide routine), it was yet considered a 
decisive advantage of the algorithm chosen that it is capable of accepting among its input, 
and of being capable to directly process, without the need for any intervening 
computational steps, the earlier mentioned tide amplitude and phase tables that are in-house 
products of NSWC. This aspect is particularly important because it is expected that, during 
the useful life of TERRA/CLLEST, the present M 2 tide tables will be augmented by 
similar tables for various other constituents of the ocean tide spectrum. Once available, a 
mere extension of the present algorithm will make possible the immediate use of the 
additional program input, without first having to obtain the additional functionalizations of 
the tidal ocean that would be needed if the second ocean tide routine had been adopted 
for inclusion in TERRA/CELEST. 
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APPENDIX A 


NOTES CONCERNING THE SECOND OCEAN TIDE ROUTINE 




NOTES CONCERNING THE SECOND OCEAN TIDE ROUTINE 


As mentioned in the INTRODUCTION, the second computer routine for the ocean 
tide is based on a functionalization of the tidal ocean surface. This functionalization 
represents the height of the tidal surface above mean sea level, h, as an expansion in 
surface harmonics. 


h = f cos (at* + x ~ 

8) = f cos 8 cos (at* + x) + f sin 8 sin (at* + x) 

(A-l) 

15 

f cos 6 = £ 

A (C,u 

cos mX + S nm sin mX)P™ (sin 0) 

(A-2) 

1 5 

f sin 6 = £ 

n = 0 

E ( C n, 

m =0 

n cos mX + S^ m sin mX)P™ (sin 6) 

(A-3) 


where 


c nm , s nm , c' nm 


X = longitude 
0 = latitude 

, S^ = expansion coefficients 

a = rate of the mean mean longitude of the moon 
X = mean mean longitude of the moon, at the beginning 
of the day UT 


Note that h is essentially a function of longitude, latitude, and time. This function is valid 
everywhere on the globe (sea and land - see the comments at the end of this appendix). 

The surface density corresponding to the tidal mass layer is 


s = ph = 


- •ja cos (at* + x) + P sin (at* + x) j 


(A-4) 


a = Gp ? cos 8 


(A-5) 


0 = Gp ? sin 8 


(A-6) 
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The gravitational potential associated with the tidal bulge is then 




. (3 dS' 


tfdS' 



—-— = T |P n (sin 6) P (sin0') 

Ir-n ~ r" +1 ' 


(A-9) 

(A-10) 


" (n - m)! , . 

+ 2 S (n + m)! P n (siin p n ( sin 9 > cos ~ * )[ (A-l 1) 


Now, evaluate 0, 


0, = GpR^ 



+ */2 
cos 6'dO 

n/2 


jfcos6_) 

(\f-r\S 


(A-12) 


Let 


sin 0 - y 

C A-13) 

sin O' - y' 

(A-14) 
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and consider the term ij from the expansion for f cos 6 , namely 


(f cos 6)^ = (C^ cos jX' + Sjj sin jX')Pj(y') 


(<^i )jj = GpR 2 ^ fdX' fdy' Pj(y') cosjX' 
o - 1 

* } E prr, [p n (y) P n (y') 

E °° (n - m)! n ) 

P » (y) 

+ GpR 2 S jj fdX' fdy' Pj(y') sin jX' 

"o 1 

i 00 R n r 

' j E |_ p „<y> p ,<y'» 

+ i?, CWPJfrtcosmlX-V)]! 


Repeated application of the orthogonality relationships for the Legendre functions, ; 
well as for the sine and cosine functions, makes possible the following equations: 

(< Mi,j=o = GpR 2 C io 27r fdy' Pj(y') V P (y) p (y') 

J i iT^o r 

= 2rrR 2 Gp 2 -^ C jo P.(y) (A-H 
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- *w’q, t ± $y’> C(y') <V 

* p^(y) fcosjx' cos m(X - X') dX'l 
0 -> 

+ 2G ' R Xt£,fS£ f4v')C(/»«y 


P n (y) J sin jX' cos m(X - X') 


JcosjX' cos m(X - X') dX' = n 5 mj cosjX 


JsinjX' cosm(X-X') dX' = tt 5 mj sinjX 


1 n = m 

if 

0 n =# m 


(0 > ) ij’ k o 2 ttCpR 2 — pi( sin Q) * (c .. cos jx + s.. sin jX) (A-22) 


and consequently. 


(4>l >« = R 2 Gp 2i+1 * (Cy cosjX+ S jj sin jX) P{(sin0) (A-23) 


Thus, the following results in a similar fashion: 


(0 2 ) (j = 2?tR 2 Gp — — * ( C ;. cos jX + S|, sin jX) P{(sin Q) (A -24) 
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Noting that the original expansion of the tidal surface is limited to n max = 15, the 
tide potential is now 


<t> = £ £ +H ,* V J 

l> = 0 M=0 


(A-25) 


27rR 2 Gp 2 r "i 

F = - -- C cos (ot* + x) + C sin (at* + x) (A-26) 

p E 2v + 1 L vtl "" J 


, 2ttR 2 Gp 2 r , -i 

H = -—— S cos (at* + x) + S sin (at* + x)J (A-27) 

2v + 1 L V>1 J 

= p E ^77 P£ (sin0) cospX (A-28) 

V v M = Fe ~^Tf p £ ( sin0 ) sin pX (A-29) 

Up to here, each of the physical quantities occurring in the equations may still be 
entered in terms of arbitrary units. However, for the purpose of the second ocean tide 
computer algorithm, it is necessary to employ compatible units: kilometers for length, 
kilograms for mass, and mean solar seconds for time. As obtained from NOAA, the 

coefficients C^, S V)j , , and are listed in meters. To accommodate the existing 
computer card deck for these data, it is necessary to multiply the right sides of Equations 

A-26 and A-27 by 10 “ 3 prior to their use in the main part of this report. 

It should finally be mentioned that the expansion for the tidal surface (Equations A-l 

through A-3) is intended to be valid for both sea and land in the sense that it is a fit to 

the NSWC/Schwiderski tidal heights and phases on the ocean and to zero tidal heights on 
land. Because of the relatively limited number of terms employed in the expansion, this fit 
is somewhat imperfect, especially near the coast lines. This imperfection is intentional. The 
reason for it is that the expansion was originally meant to be the basis for finding a tide 
potential for satellite work rather than a continuous, highly accurate description of the 

tidal surface. Satellite orbits appear to be much more sensitive to the low-order harmonics 
in the tidal disturbing potential than to the higher ones. Based on information from an 
outside source, we in fact anticipate that certain of the second- and fourth-order harmonics 
will be the main contributors to the orbital perturbation. Thus, the emphasis on the 
low-order terms apparent in the expansion for the tidal surface is rather an advantage than 
a defect as far as our own satellite work is concerned. 
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APPENDIX B 

A COMPARISON OF THE PERTURBING ACCELERATION DUE TO THE 
SOLID EARTH TIDE WITH THAT CAUSED BY THE OCEAN TIDE 
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A COMPARISON OF THE PERTURBING ACCELERATION DUE TO THE 
SOLID EARTH TIDE WITH THAT CAUSED BY THE OCEAN TIDE 


Preliminary runs of the first computer routine for the ocean tide and of the routine 
for the solid earth tide, all of which used complete sets of input data (as opposed to the 
restricted input data of the numerical trial cases listed in the main body of this reporti. 
have resulted in the following perturbing accelerations. 

The satellite coordinates and the perturbing acceleration components listed are referred 
to the inertial frame corresponding to the mean equinox and mean equator of J = 1977. 
n = 202, and t* = 0 sec. 


Time t 


J = 1977 
n = 202 
t* = 50000 sec 

Satellite coordinates at t 

x, = 3151.52923 km 
x 2 = 5458.60875 
x 3 = 3639.07250 


Perturbing acceleration due to ocean tide at t 


T x = + 0.594629E - 11 km sec 2 
T y = + 0.134622E - 10 
T z = - 0.258193E - 11 


Perturbing acceleration due to solid earth tide at t 


T x = + 0.149368E - 09 km sec -2 
T y = + 0.125475E - 10 
T z = - 0.184573E - 10 

T x = - 0.493523E - 10 
T y = +0.381372E- 10 
T z = + 0.713397E - 12 


Nominal values for 
Love numbers 
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T x = + 0.164304H - 09 
T y = + 0.013802E - 09 
T z = -0.020303E - 09 

T x = -0.054288E - 09 
T* = + 0.041951E - 09 
T z = + 0.000785E- 09 

T x = + 0.179241E - 09 
T y = + 0.150570E - 10 
T z = - 0.221487E - 10 

T x = -0.592228E - 10 
T y = + 0.457647E - 10 
T z = + 0.856076E - 12 


Each Love number 
increased by 10 percent 


due to moon 


Each Love number 
increased by 20 percent 


The magnitudes of the perturbing accelerations are then 

T(Ocean tide) = T QCEAN = 1.49E- 11 km sec -2 

T(Solid earth tide due to moon, for nominal Love Numbers) 

= T SOL)d = 1.51026E- 10 km sec' 2 

T(Solid earth tide due to moon, for Love numbers increased by 10 percent) 
= T solid + 10 = 1.66048E- 10 km sec -2 

T(Solid earth tide due to moon, for Love numbers increased by 20 percent) 
= T so l i d + 2 o = 1-8123IE- 10 km sec' 2 

Note that 


^SOLID+I0 “ ^SOLID * ^OCEAN (B-l) 

or 

^SOLID + 10 '^SOLID I 55 I^OCEAnI (®*2) 
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This indicates that, for the particular combination of time and satellite position 
chosen, if one increases the Love numbers by ten percent, the perturbing acceleration due 
to the lunar component of the solid earth tide increases by an amount roughly equal to 
the perturbing acceleration due to the ocean tide. For the particular example chosen, the 

two equations are accurate to better than two percent, which is suspected to be a 

coincidence. But we are certain that this example reveals the order of magnitude of a 

relationship existing between the gravitational effects of the ocean tide and the solid earth 

tide. 


73 



DISTRIBUTION 


Defense Mapping Agency, Headquarters 
Naval Observatory Building 56 
Washington, DC 20360 

Defense Mapping Agency 
Aerospace Center 
St. Louis, MO 63118 
ATTN: Richard L. Ealum 

Defense Mapping Agency 
Hydrographic Center 
Washington, DC 20390 

Defense Mapping Agency 
Topographic Center 
6500 Brooks Lane 
Washington, DC 20315 

Chief of Naval Research 
Department of the Navy 
Arlington, VA 22217 

Department of the Air Force 
HQ Space and Missile Test Center (AFSC) 
Vandenberg Air Force Base, CA '93437 
ATTN: Sven C. Ernberg 

Naval Research Laboratory 
Washington, DC 20360 
ATTN: Code 4130 

U.S. Navy Astronautics Group 
Point Mugu, CA 93041 
ATTN: Code AG33 
Code SPMOO 
Mr. Campbell 

National Oceanic and Atmospheric Administration 
National Ocean Survey 
Rockville, MD 20850 



National Aeronautics and Space Administration 
Washington, DC 20546 
ATTN: KSS-10 


Aerospace Corporation 
2350 East El Segundo Boulevard 
El Segundo, CA 90245 
ATTN: Library 

Applied Physics Laboratory 
The Johns Hopkins University 
8621 Georgia Avenue 
Silver Spring, MD 20910 

Library of Congress 
Washington, DC 20540 

ATTN: Gift and Exchange Division (4) 

Defense Documentation Center 
Cameron Station 

Alexandria, VA 22314 (12) 

Local: 

E41 

E416 (Cannon) 

KOI 

K05 (Schwiderski) 

K10 

K11 (Groeger) 

K50 
K70 
V 

X210 

X2101 (GIDEP) 

X211 


(50) 

(30) 

(5) 

(5) 

( 2 ) 

( 2 ) 



